https://rspatial.org/raster/rs/2-exploration.html

1 Data

1.2 Unzip data

unzip('data/rsdata.zip', exdir='data')
## Warning in unzip("data/rsdata.zip", exdir = "data"): error 1 in extracting from
## zip file

2 Image properties

“Create RasterLayer objects for single Landsat layers (bands)”

library(raster)
## Loading required package: sp
# Blue
b2 <- raster('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- raster('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- raster('data/rs/LC08_044034_20170614_B4.tif')
# Near Infrared (NIR)
b5 <- raster('data/rs/LC08_044034_20170614_B5.tif')

2.1 Check variables

see the spatial resolution, extent, number of layers, coordinate reference system, etc

b2
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : LC08_044034_20170614_B2.tif 
## names      : LC08_044034_20170614_B2 
## values     : 0.0748399, 0.7177562  (min, max)
b3
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : LC08_044034_20170614_B3.tif 
## names      : LC08_044034_20170614_B3 
## values     : 0.04259216, 0.6924697  (min, max)
b4
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : LC08_044034_20170614_B4.tif 
## names      : LC08_044034_20170614_B4 
## values     : 0.02084067, 0.7861769  (min, max)
b5
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : LC08_044034_20170614_B5.tif 
## names      : LC08_044034_20170614_B5 
## values     : 0.0008457669, 1.012432  (min, max)

2.1.1 Check CRS (coordinate reference system)

# coordinate reference system (CRS)
crs(b2)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b3)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b4)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b5)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
# Number of cells, rows, columns
ncell(b2)
## [1] 1863765
dim(b2)
## [1] 1245 1497    1
# spatial resolution
res(b2)
## [1] 30 30
# Number of bands
nlayers(b2)
## [1] 1
# Do the bands have the same extent, number of rows and columns, projection, resolution, and origin
compareRaster(b2,b3)
## [1] TRUE

3 Create raster stack object of multiple layers from single raster band objects

s <- stack(b5, b4, b3)
s
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B5, LC08_044034_20170614_B4, LC08_044034_20170614_B3 
## min values :            0.0008457669,            0.0208406653,            0.0425921641 
## max values :               1.0124315,               0.7861769,               0.6924697
plot(s)

plot(b5)

3.1 Raster stacks created using filenames

3.1.1 List filenames

# create list of raster objects
filenames <- paste0("data/rs/LC08_044034_20170614_B", 1:11, ".tif")
filenames
##  [1] "data/rs/LC08_044034_20170614_B1.tif" 
##  [2] "data/rs/LC08_044034_20170614_B2.tif" 
##  [3] "data/rs/LC08_044034_20170614_B3.tif" 
##  [4] "data/rs/LC08_044034_20170614_B4.tif" 
##  [5] "data/rs/LC08_044034_20170614_B5.tif" 
##  [6] "data/rs/LC08_044034_20170614_B6.tif" 
##  [7] "data/rs/LC08_044034_20170614_B7.tif" 
##  [8] "data/rs/LC08_044034_20170614_B8.tif" 
##  [9] "data/rs/LC08_044034_20170614_B9.tif" 
## [10] "data/rs/LC08_044034_20170614_B10.tif"
## [11] "data/rs/LC08_044034_20170614_B11.tif"

3.1.2 Create an 11 layer raster stack from all 11 filenames

“The layers represent reflection intensity in the following wavelengths: Ultra Blue, Blue, Green, Red, Near Infrared (NIR), Shortwave Infrared (SWIR) 1, Shortwave Infrared (SWIR) 2, Panchromatic, Cirrus, Thermal Infrared (TIRS) 1, Thermal Infrared (TIRS) 2.
We won’t use the last four layers and you will see how to remove those in following sections.”

landsat <- stack(filenames)
landsat
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 11  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11 
## min values :            9.641791e-02,            7.483990e-02,            4.259216e-02,            2.084067e-02,            8.457669e-04,           -7.872183e-03,           -5.052945e-03,            3.931751e-02,           -4.337332e-04,             2.897978e+02,             2.885000e+02 
## max values :              0.73462820,              0.71775615,              0.69246972,              0.78617686,              1.01243150,              1.04320455,              1.11793602,              0.82673049,              0.03547901,             322.43139648,             317.99530029

4 Maps

4.1 Grey-scale Single band and composite maps

plot individual layers in a 2x2 grid
legend values between 0-1
“different surface features reflect the incident solar radiation differently.
Each layer represent how much incident solar radiation is reflected for a particular wavelength range.
For example, vegetation reflects more energy in NIR than other wavelengths and thus appears brighter.
In contrast, water absorbs most of the energy in the NIR wavelength and it appears dark.”

par(mfrow = c(2, 2))
plot(b2, main = "Blue", col = gray(0:100 / 100))
plot(b3, main = "Green", col = gray(0:100 / 100))
plot(b4, main = "Red", col = gray(0:100 / 100))
plot(b5, main = "NIR", col = gray(0:100 / 100))

4.2 True / natural colour image

“(vegetation in green, water blue etc), we need bands in the red, green and blue regions.
For this Landsat image, band 4 (red), 3 (green), and 2 (blue) can be used.”

par(oma = c(2, 0, 1, 4)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(0.8, 3.1, .55, 0.35))

landsatRGB <- stack(b4, b3, b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")

## False colour image

par(oma = c(2, 0, 1, 4)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(0.8, 3.1, .55, 0.35))

landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")

4.3 Compare True and False Colour Images

par(mfrow = c(1, 2))
par(oma = c(0, 0.5, 0, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "s") # "s" square plotting region vs m - maximal


plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")

landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")

5 Arguements to modify image

help(plotRGB)

6 Subset (indexing) and rename bands in landsat object of all layers

# select first 3 bands
landsat
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 11  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11 
## min values :            9.641791e-02,            7.483990e-02,            4.259216e-02,            2.084067e-02,            8.457669e-04,           -7.872183e-03,           -5.052945e-03,            3.931751e-02,           -4.337332e-04,             2.897978e+02,             2.885000e+02 
## max values :              0.73462820,              0.71775615,              0.69246972,              0.78617686,              1.01243150,              1.04320455,              1.11793602,              0.82673049,              0.03547901,             322.43139648,             317.99530029
landsatsub1 <- subset(landsat, 4:2)
landsatsub1
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B4, LC08_044034_20170614_B3, LC08_044034_20170614_B2 
## min values :              0.02084067,              0.04259216,              0.07483990 
## max values :               0.7861769,               0.6924697,               0.7177562
# same as
landsatsub2 <- landsat[[4:2]]
landsatsub2
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B4, LC08_044034_20170614_B3, LC08_044034_20170614_B2 
## min values :              0.02084067,              0.04259216,              0.07483990 
## max values :               0.7861769,               0.6924697,               0.7177562
# Number of bands in the original and new data
nlayers(landsat)
## [1] 11
nlayers(landsatsub1)
## [1] 3
nlayers(landsatsub2)
## [1] 3

6.1 Remove last 4 bands that won’t be used

landsat <- subset(landsat, 1:7)
landsat
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7 
## min values :            0.0964179114,            0.0748399049,            0.0425921641,            0.0208406653,            0.0008457669,           -0.0078721829,           -0.0050529451 
## max values :               0.7346282,               0.7177562,               0.6924697,               0.7861769,               1.0124315,               1.0432045,               1.1179360

6.2 Rename Band names

names(landsat)
## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
names(landsat)
## [1] "ultra.blue" "blue"       "green"      "red"        "NIR"       
## [6] "SWIR1"      "SWIR2"

6.3 Spatial subset or crop

“used to limit analysis to a geographic subset of the image”

extent(landsat)
## class      : Extent 
## xmin       : 594090 
## xmax       : 639000 
## ymin       : 4190190 
## ymax       : 4227540
e <- extent(624387, 635752, 4200047, 4210939)
e
## class      : Extent 
## xmin       : 624387 
## xmax       : 635752 
## ymin       : 4200047 
## ymax       : 4210939
# crop landsat by the extent
landsatcrop <- crop(landsat, e)
landsatcrop
## class      : RasterBrick 
## dimensions : 363, 379, 137577, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 624390, 635760, 4200060, 4210950  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      :  ultra.blue,        blue,       green,         red,         NIR,       SWIR1,       SWIR2 
## min values : 0.101796150, 0.075989284, 0.044045158, 0.030556191, 0.007243267, 0.001865030, 0.003144530 
## max values :   0.4548080,   0.4746728,   0.5034941,   0.5592065,   0.7013395,   0.9038041,   0.9750224

7 Compare landsat True and False Colour Images

par(mfrow = c(1, 2))

landsatsub3 <- subset(landsat, 4:2)
landsatsub4 <- subset(landsat, 5:3)

plotRGB(landsatsub3, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")

plotRGB(landsatsub4, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")

8 Compare cropped landsat True and False Colour Images

8.1 Cropped Landsat True Colour Image

par(ps = 9) # point text size
par(oma = c(0, 1.5, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(1.2, 4.2, .55, 0.35))
par(mar = c(5, 5, 4, 2))
# ?par # help for par() function

landsatsub5 <- subset(landsatcrop, 4:2)
plotRGB(landsatsub5, axes = TRUE, xlab = "UTM E/W", ylab = "UTM N/S", stretch = "lin", main = "Landsat\nTrue Colour Composite")

8.2 Cropped Landsat False Colour Image

par(ps = 9) # point text size
par(oma = c(0, 1.5, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(1.2, 4.2, .55, 0.35))
par(mar = c(5, 5, 4, 2))
# ?par # help for par() function

# title(sub = "Comparing True and False Colour Landsat Composites") # only works for area of one graph
# xlab = "UTM E/W", ylab = "UTM N/S",  takes too much space
landsatsub6 <- subset(landsatcrop, 5:3)
plotRGB(landsatsub6, axes = TRUE, stretch = "lin", main = "Landsat\nFalse Colour Composite")

8.3 Compare cropped landsat True and False Colour Images

par(mfrow = c(1, 2))
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
# ?par # help for par() function

landsatsub5 <- subset(landsatcrop, 4:2)
plotRGB(landsatsub5, axes = TRUE, xlab = "UTM E/W", ylab = "UTM N/S", stretch = "lin", main = "Landsat\nTrue Colour Composite")

# title(sub = "Comparing True and False Colour Landsat Composites") # only works for area of one graph
landsatsub6 <- subset(landsatcrop, 5:3)
plotRGB(landsatsub6, axes = TRUE, xlab = "UTM E/W", ylab = "UTM N/S", stretch = "lin", main = "Landsat\nFalse Colour Composite")

9 Saving results to disk

9.1 Save as jpg

jpeg(file="landsatRGB-TCC.jpeg")
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
dev.off()
## quartz_off_screen 
##                 2
jpeg(file="landsatRGB-FCC.jpeg")
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
dev.off()
## quartz_off_screen 
##                 2

9.2 Save plot as tiff

tiff(file="landsatRGB-TCC.tif")
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
dev.off()
## quartz_off_screen 
##                 2
tiff(file="landsatRGB-FCC.tif")
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
dev.off()
## quartz_off_screen 
##                 2

9.3 Save writeRaster as tiff

rf <- writeRaster(landsatcrop, filename = "cropped-landsat.tif", overwrite=TRUE)
rf <- writeRaster(landsatsub5, filename = "cropped-landsat5.tif", overwrite=TRUE)
rf <- writeRaster(s, filename = "cropped-landsat-TCC.tif", overwrite=TRUE)

# S4 method for RasterStackBrick,character

rf <- writeRaster(landsatFCC, filename=file.path("landsatFCC.tif"), format="GTiff", overwrite=TRUE)
rf <- writeRaster(landsatRGB, filename=file.path("landsatTCC.tif"), format="GTiff", overwrite=TRUE)

class(landsatcrop)
## [1] "RasterBrick"
## attr(,"package")
## [1] "raster"
class(landsat)
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"

10 Relationships between raster layer Bands

10.1 Plot of UV/B reflection correlation

There is a high correlation (99%) between the ultra-blue and blue wavelengths.
so not much info will be lost if sing only one of these bands

pairs(landsatcrop[[1:2]], main = "Correlation between Ultra-blue versus Blue")

10.2 Plot of R/NIR wavelength reflection Correlation

A triangle distribution. Vegetation reflects higher in NIR than red wavelengths
“and creates the upper corner close to NIR (y) axis.
Water absorbs energy from all the bands and occupies the location close to origin.
The furthest corner is created due to highly reflecting surface features like bright soil or concrete.”

pairs(landsatcrop[[4:5]], main = "Red vs NIR")

11 Land Cover

11.1 Extract pixel values from land cover info polygons

# load the polygons with land use land cover information
samp <- readRDS('data/rs/samples.rds')

# generate 300 point samples from the polygons
ptsamp <- spsample(samp, 300, type='regular')
## Warning in proj4string(obj): CRS object has comment, which is lost in output
# add the land cover class to the points
ptsamp$class <- over(ptsamp, samp)$class

# extract values with points
df <- extract(landsat, ptsamp)

# To see some of the reflectance values
head(df)
##      ultra.blue      blue      green        red       NIR     SWIR1     SWIR2
## [1,]  0.1294681 0.1235693 0.12922950 0.17581198 0.3241255 0.3347085 0.1969346
## [2,]  0.1422847 0.1373185 0.14692563 0.18604797 0.2774130 0.3170558 0.2110958
## [3,]  0.1327861 0.1148514 0.09830463 0.09570226 0.1651856 0.2010984 0.1640145
## [4,]  0.1344559 0.1170851 0.10006124 0.09889016 0.1687422 0.2109006 0.1755517
## [5,]  0.1311379 0.1291861 0.14445338 0.19541651 0.3324314 0.3122631 0.1759638
## [6,]  0.1384245 0.1338053 0.14467025 0.19270571 0.3112655 0.3450096 0.2207463

11.2 Spectral profiles of land cover

ms <- aggregate(df, list(ptsamp$class), mean)

# instead of the first column, we use row names
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms
##          ultra.blue       blue      green        red        NIR      SWIR1
## built     0.1749001 0.16639696 0.16668714 0.17878095 0.22796993 0.23626448
## cropland  0.1118587 0.08961921 0.08389110 0.05294961 0.46342623 0.15512600
## fallow    0.1321714 0.11686468 0.10508105 0.11573817 0.18102815 0.23419503
## open      0.1385035 0.13693335 0.15151051 0.20491259 0.34288999 0.35559024
## water     0.1337126 0.11681711 0.09981322 0.07909337 0.04886901 0.03346269
##               SWIR2
## built    0.19687054
## cropland 0.06828988
## fallow   0.19450095
## open     0.21181017
## water    0.02714454

11.2.1 Plot mean spectral Profile of features

"The spectral profile plots reflectance of different earth’s surface features
‘Water’ has low reflection in all wavelengths
‘built’, ‘fallow’ and ‘open’ have have high reflectance in longer wavelengths

# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')

#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)

# First create an empty plot
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")

# add the different classes
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}

# Title
title(main="Spectral Profile from Landsat", font.main = 2)

# Legend
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")

12 Basic Math with raster

performed on each pixel / grid cell

12.1 Combine bands

12.1.1 Data (same as above)

raslist <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <-stack(raslist)
landsatRGB <- landsat[[c(4, 3, 2)]]
landsatFCC <- landsat[[c(5, 4, 3)]]

12.1.2 Function Vegetation Index

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}

12.1.2.1 NDVI (Normalized Difference Vegetation Index)

to view greenness
Landsat bands NIR = 5, red = 4 #### NDVI one way

ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Landsat NDVI", sub = "Level of greenness")

12.1.2.2 NDVI another function

vi2 <- function(x, y){
  (x - y) / (x + y)
}
ndvi2 <- overlay(landsat[[5]], landsat[[4]], fun = vi2)
plot(ndvi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Landsat NDVI", sub = "Level of greenness")

12.1.3 Adapt code for other indices

see spectral profile for bands with max / min reflectance

12.1.3.1 Function for 2 band index

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}

12.1.3.2 Water body index NDWI (Normalized difference water index)

Landsat bands: green = 3 and NIR = 5

ndwi <- vi(landsat, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Water Body Index: Landsat Bands green (3) and NIR (5)", sub = "Monitor water in water bodies")

12.1.3.3 Water in leaf index NDWI (Normalized difference water index)

Landsat bands: NIR = 5 and SWIR = 7
(Note: very similar to Built Environment Index 2 Bands: NIR = 5 and green = 3)

ndwi2 <- vi(landsat, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Leaf Water Index: Landsat Bands NIR (5) and SWIR (7)", sub = "Monitor change in leaf water content")

12.1.3.4 Another try at Water in leaf index NDWI (Normalized difference water index)

Landsat bands: NIR = 5 and SWIR = 6
this is not the band to identify water in leaves

ndwi2 <- vi(landsat, 5, 6)
plot(ndwi2, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Leaf Water Index: Landsat Bands NIR (5) and SWIR (6)", sub = "Monitor change in leaf water content")

12.1.4 Adapt code for other indices

see spectral profile for bands with max / min reflectance
e.g. Identify built areas NDBI (Normalized Difference Built-up Index)
Landsat bands: 4 = red, 2 = green, 7 = SWIR
https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm

12.1.4.1 Function for Built Environment Index

# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
  bk <- img[[k]]
  bi <- img[[i]]
  bp <- img[[p]]
  bi <- (bk - bi) * bp / (bk + bi) * bp
  return(bi)
}

12.1.4.2 Built Environment Index 3 Bands

Landsat Bands: 4 = red, 2 = blue, 7 = SWIR
https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm

ndbi <- bi(landsat, 4, 2, 7)
plot(ndbi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands red (4), blue (2), and SWIR (7)", sub = "Monitor change in built environment (4-2)*7/(4+2)*7")

12.1.4.3 Built Environment Index 2 Bands SWIR and NIR

Landsat bands: 6 = SWIR, 5 = NIR
https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm

# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  
  bi <- (bk - bi) / (bk + bi) 
  return(bi)
}

ndbi <- bi(landsat, 6, 5)
plot(ndbi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands SWIR (6) and NIR (5)", sub = "Monitor change in built environment (6-5)/(6+5)")

12.1.4.4 Built Environment Index 2 Bands red and NIR

Landsat bands: 4 = red, 5 = NIR
** Don’t Use (Note: too similar to Built Environment Index 2 Bands: NIR = 5 and green = 3) https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm

bei <- vi(landsat, 4, 5)
plot(bei, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands NIR (5) and green (3)", sub = "Monitor change in built environment")

13 Comparison Plots

13.1 Compare Land Use Change Plots

par(mfrow = c(2, 2)) # rows, columns
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(0.35, 0.35, 0.35, 0.35))
# ?par # help for par() function

###############
# NDVI plot
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat NDVI Greeness Index")

###############
# NDWI plot leaf water
ndwi2 <- vi(landsat, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), main = "Leaf Water Index: Landsat Bands NIR (5) and SWIR (7)")

###############
# NDWI plot water bodies
ndwi <- vi(landsat, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), main = "Water Body Index: Landsat Bands green (3) and NIR (5)")

###############
# Built environment index
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
  bk <- img[[k]]
  bi <- img[[i]]
  bp <- img[[p]]
  bi <- (bk - bi) * bp / (bk + bi) * bp
  return(bi)
}

bsi <- bi(landsat, 4, 2, 7)
plot(bsi, col = rev(terrain.colors(10)), main = "Built Index: Landsat Bands Red, Green and SWIR [(4-2)*7/(4+2)*7]")

13.2 Compare cropped image Land Use Change Plots

par(mfrow = c(2, 2)) # rows, columns
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(0.35, 0.35, 0.35, 0.35))
# ?par # help for par() function

###############
# NDVI plot
ndvi <- vi(landsatcrop, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat NDVI Greeness Index")

###############
# NDWI plot leaf water
ndwi2 <- vi(landsatcrop, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), main = "Leaf Water Index: Landsat Bands NIR (5) and SWIR (7)")

###############
# NDWI plot water bodies
ndwi <- vi(landsatcrop, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), main = "Water Body Index: Landsat Bands green (3) and NIR (5)")
###############
# Built environment index
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
  bk <- img[[k]]
  bi <- img[[i]]
  bp <- img[[p]]
  bi <- (bk - bi) * bp / (bk + bi) * bp
  return(bi)
}

bsi <- bi(landsatcrop, 4, 2, 7)
plot(bsi, col = rev(terrain.colors(10)), main = "Built Index: Landsat Bands Red, Green and SWIR [(4-2)*7/(4+2)*7]")

13.3 Spectral Plot

###############
# Plot Spectra
# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')

#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)

# First create an empty plot
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")

# add the different classes
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}

# Title
title(main="Spectral Profile from Landsat", font.main = 2)

# Legend
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")